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The development of computational fluid dynamics over the last few decades has yielded enormous successes 
and capabilities being routinely employed today; however there remain some open problems to be properly 
resolved-some are fundamental in nature and some resolvable by operational changes. These two categories 
are distinguished and broadly explored previosuly. * 1 One, that belongs to the former, is the so-called overheat- 
ing problem, especially in rarefying flow. This problem up to date still dogs every method known to the author; 
a solution to it remains elusive. The study in this paper concludes that: (1) the entropy increase is quantita- 
tively linked to the increase in the temperature increase, (2) it is argued that the overheating is inevitable 
in the current shock capturing or traditional finite difference framework, and (3) a simple hybrid method is 
proposed that removes the overheating problem in the rarefying problems, but also retains the property of 
accurate shock capturing. This remedy (enhancement of current numerical methods) can be included easily in 
the present Eulerian codes. 


I. Introduction 

T he overheating problem remains one of the most prominent open CFD problems; it is typically found in two sce- 
nariosone related to colliding flows and the other to receding flows. The former, resulting either from reflection 
of a flow from a wall or equivalently from two colliding (of equal or different strength) and manifested as a spike 
in temperature near the wall, is a well-known one since it was first recognized by von Neumann. 2 This overheating 
problem still vexes CFD researchers today (e.g., books by Toro 3 and LeVeque 4 ). The latter, while less known, occurs 
during rarefaction. For example, when two streams receding in relative motion from each other, a rarefying region 
between them is being generated, where all thermodynamic propertiespressure, density, and temperature arc decreas- 
ing. However, numerical solutions consistently show an increased temperature even though the other two variables 
behave correctly, as shown in Fig. 1 in which two typical solutions are displayedone by a sharp-shock method (such as 
the Godunov method) and the other by a diffusive-shock method (e.g. HLLE method).Both cases, although resulting 
from different scenarios, share a common numerical feature of elevated temperature: it appears in both the Eulerian 
and the Lagrangian approaches.? But there is a fundamental difference between them: the receding flow should exert 
no generation of entropy, while the colliding flow naturally is expected to see an increased entropy at the instance of 
impact. 


Our study of this long-standing problem is part of our broader interest in exploring open problems in shock captur- 
ing, specifically those typically adopted in research or production environments. This arguably limits to the Eulerian 
formulation and to a large extent, to finite volume discretization. The exploratory study, first reported in the previous 
AIAA CFD Conference, 1 identified that the overheating in rarefying flow is most challenging, defies physics (in the 
form of temperature increase), and remains an unsolved mystery. A follow-up study 5 further attempted to find the 
cause of the overheating, discovering the violation of entropy constancy in this rarefying flow and its close connec- 
tion to overheating. Our interest is intensified by the desire to understand the numerical mechanism responsible for 
the overheating on the one hand, but also by the challenge of mitigating and correcting this error. The current paper 
reports further progress in the following aspects: (1) the entropy increase is quantitatively linked to the increase in 
the temperature increase, (2) it is argued that the overheating is inevitable in the current shock capturing or traditional 
finite difference framework, and (3) a simple hybrid method is proposed that removes the overheating problem in the 
rarefying problems, but also retains the property of accurate shock capturing. This remedy (enhancement of current 
numerical methods) can be easily included in the present Eulerian codes. 

^Senior Technologist, Aeropropulsion Division, MS 5-11, 21000 Brookpark Road, Cleveland, OH 44135. Associate Fellow. 
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Figure 1. Receding flows computed on a grid of 100 cells, with the initial condition, (p,p, u)l = (1.0, 0.4, —2.0) and (p, p, u)h = (1.0, 0.4, 2.0). 


II. Problem statement: Receding flows 

In the discussion of the overheating problem and its remedy, we shall exclusively focus on the receding problem. 
The simplest generic form for the receding problem is a Riemann problem with the initial condition: 

u L <0 >Pl,Pl and u R > 0,p R} p R (1) 

where the “L” and “R” states may be of same or different strength. The strength of rarefaction is defined by the relativel 
Mach number | Mr - Ml |; different combinations of thermodynamic properties will of course result in different flow 
configurations. 

As stated above, the computation will be carried out in the finite-volume Eulerian setting. The exact solutions 
as they are available will be used as a reference for comparison with calculated solutions. We will use the Godunov 
result 3 both to establish the baseline solution for comparison with other methods, and also to gauge its own ability 
to solve the considered problems against the exact solution. The explicit Euler and two-step Runge-Kutta methods 
are used for the first- and second-order time integrations respectively. To evaluate numerical methods, we ensure that 
solutions are physically relevant, numerically stable and mesh-size converged. Hence, it may be necessary with some 
methods to reduce the time step (CFL) further than others. In other words, a non-converged numerical solution will not 
be considered as a valid solution of the problem studied. Subsequently, results obtained by the AUSM + -up method 6 
will be presented for comparison purposes. 
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(a) M=0.8, 100 cells. 


(c) M=2.5, 100 cells. 


(d) M=2.5, 10,000 cells. 


(b) M=0.8, 10,000 cells. 


Figure 2. Solutions of receding flows moving at M=0.8 and 2.5 by the Godunov method on 100 and 10,000 cells respectively. 


In Fig. 2, the overheating problem exists in both low and high speeds and the Godunov solutions reveal: (1) the 
overheating error from the true solution remains unchanged after increasing the number of cells by 100 times and (2) 
the error increases with the magnitude of the separation velocity (the strength of rarefaction). A diffusive method, like 
HLLE, clearly spreads the error in spatial direction. 

a The Godunov method is used as a reference as it has been regarded as the gold standard. 
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In our previous study, 1 results from other methods (such as van Leer, HLLE, SLAU2, AUSM^-up) are evaluated, 
with extensive studies on order of accuracy and grid convergence. However, the cause of the overheating problem 
remained unknown and no clear clue as to how it can be overcome emerged. Some representative results are given 
here to provide a frame of reference, they are summarized below. 

II.A. Effect of spatial discretization accuracy and limiters 

Despite the fact that increasing the number of grid points does not lead a way out of the overheating problem, it 
would be interesting to see if higher order discretization might be enlightening, since the solution indeed converges 
in the major part of the domain in Fig. 2. In Fig. 3 we compare the second-order accurate results with the first-order 
one, by using primitive and characteristic variables for interpolation together with the diffusive and sharp-minmod and 
superbee limiters. We observe: (1) the stable time step has to be significantly reduced, (2) the use of primitive variables 
has limited success in reducing the overheating, (3) a significant improvement arises from the use of characteristic 
variables, and (4) the minmod and superbee limiters have a dramatic and opposite influence on the entropy distribution 
and the temperature distribution in the middle low speed region. Notice that the superbee limiter causes a decrease in 
entropy, thereby producing an overcooling instead! In fact the solution was so overcooled during the calculation that a 
negative density occurred and the solution was artificially kept alive by resetting the density to a small positive value. b 
Thus, we reiterate: 

Observation 1: The “overheating” by the first-order method is not reduced by refining the size of grid or time step . 
A significant reduction in overheating can be achieved when using the second-order interpolation ; but the opposite, 
overcooling in this case, can be produced by combining the superbee limiter with the characteristic variables. Hence, 
the problem remains resolved . 
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Figure 3. solutions computed by the Godunov method with the first and second order accuracy together with various limiter functions and 

interpolated variables with 100 cells. 


II.B, Effect of numerical flux function 

We extensively examined how methods other than the Godunov behaved and presented the findings in previous pa- 
pers. 1,5 Here we will merely show the representative results by the AUSM+-up method. Figure 4 shows that the 
accuracy of these results is comparable to their Godunov counterparts. The peak entropy values by the first order and 
second-order methods are nearly identical to each other respectively, with insignificant differences. 

b This result should have been discounted and it is included to highlight the need for investigating closely the computed result in every aspect, 
even though the temperature appears to be in a correct trend. 
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Figure 4. M^=-M«=-2.5; solutions computed by the AUSM + -up method with the first and second order accuracy together with various limiter functions and 
interpolated variables with 100 cells. 


These results give the first hint of a close connection between the entropy production and overheating. It also 
implies that numerical diffusion, introduced by the first and second order discretization, plays a significant role in 
reducing, if not correcting, the overheating anomaly. Both the results from using the characteristic variables with the 
minmod and superbee limiters shown in Figs. 3 and 4 give the first hint about the role of entropy. 

Observation 2: The entropy is a good indicator for diagnosis because its value from computation varies proportionally 
with the temperature , while theoretically the entropy remains constant and the temperature drops during receding. 

In a follow-up study, 5 we identified the link between the overheating and the entropy generation (an anomaly 
for this receding problem) and proved the necessity of preserving the entropy constancy in order to obtain a correct 
solution. This finding provides the basis for studying the entropy distribution. It was anticipated the key for unlocking 
the difficulty would lie in the ability to eliminate the entropy generation. Hence, in the present study, we focus on the 
entropy function. 


III. Entropy function: Why is it useful, but not used? 


To find how the temperature rise begins, we need to see the evolution of the flowfield. Figure 5 displays the 
temperature distribution at several time slices, in which the early times are plotted on the left figure in a blow-up scale 
and the later times are on the right. We see that the temperature just within the first time step (shown in blue) is already 
higher than it should be. We also remark that for this receding problem while the pressure and density are vanishing, 
the temperature reaches a non-zero constant value. Figure 6 shows that the entropy is also generated in the first time 
step and continues to grow spatially and temporally. More strikingly, Fig. 7 shows the peak entropy attains a saturated 
value, which remains unchanged even the grid size is refined, suggesting that the root of overheating lies out reach of 
simple fixes by numerical discretization alone. 
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Figure 5. Evolution of internal energy, plotted at various time slices. 


4 of 13 


American Institute of Aeronautics and Astronautics 




X 


X 


o 

</) 

(/) 


Figure 6. Evolution of entropy at various time slices. 
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Figure 7. Saturated entropy on different grids. 


The fundamental thermodynamic equation applied to an ideal gas can be written as 

Tds = de- -j^dp = e(^~ -7^) = ed(log^) = C V T d(log ~). 


( 2 ) 


Hence, the last equality allows us to define the entropy variable for an ideal gas as (omitting the proportional constant 
C v without loss of generality) as 

s = 1 °gp7 (3) 

Rewriting the energy balance equation with the substitution of continuity and momentum equations, we find 

De p Dp 

Dt _ p2Of =0 ' (4) 

Comparing Eqs. 2 and 4, we get the entropy transport equation (under the assumption of smooth flow), 

Ds ds ds 

m~di + Ui d^i~ 0 ' (5) 

Combining the first law of thermodynamics with Eq. 2 and expressing the work dW = -pd(l/p), we have 7 

dq = Tds. (6) 
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Assuming the heat input to the ’’numerical” system dq is generated via numerical diffusion, 

dq = KdT (7) 

here k is interpreted as the numerical heat conductivity, assumed to be constant in our analysis. Substituting into Eq. 
6 and integrating, we get 

As = s(x (2) )-s(x (1) ) = /dog^^ (8) 

If the heat input that is to contribute the overheating does so via heat diffusivity, then the numerical entropy should 
follow the above rule in the numerical solution: (1) the relationship between the entropy and temperature is one to 
one; the entropy increases as the temperature increases and (2) these two quantities follow the above logarithmic 
relationship. In Fig. 9 we plot the variation of these two quantities between two neighboring points for the entire 
computation domain, it reveals that the relationship arrived at in Eq. 8 holds, confirming that entropy increases with 
temperature and the temperature increase is activated through numerical diffusion. In the region where the numerical 
rarefaction wave has not reached, there are no changes in s and T and they are clustered at the center of the plot; the 
points further away from the center correspond to the region of large changes in entropy and temperature, namely the 
center of the receding flow. 

On the other hand, using the entropy equation Eq. 5, in place of the conventional energy equation, completely 
eradicates the overheating and yields a correct and grid-converged solution, as shown in Fig. 8. Hence, one is led to 
confirm directly if by enforcing the entropy conservation the overheating problem is solved. 










(a) 100 cells, 0(h). (b) 400 cells, 0(h). 

Figure 8. Solution of the receding flow by including the entropy transport equation in place of the energy equation, on 100-cell and 400-cell grids, showing 
elimination of overheating and grid convergence. 


Equation 5 states that the entropy remains constant following the particle with velocity u\ it describes the con- 
vection of an initial entropy profile by u(x,t). However, nowhere does the entropy function appear explicitly in the 
continuity and momentum equations, while the pressure can be directly determined from s by knowing the density 
from the continuity equation. The entropy convection gives no mechanism for describing an entropy-increasing event, 
such as shock waves. That is, the entropy function remains constant across the shock wave in a shock tube problem as 
will be evident later, thus rendering the entropy equation unusable for general problems. 
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Figure 9. Temperature variation A x log T vs. entropy variation A x s. 


Moreover, the second law of thermodynamics only demands a nondecreasing change in entropy ( ds > 0), the 
strength of this inequality is clearly problem-dependent and the upper bound is not known a priori. Since a pre- 
cise mathematical expression for describing the convection and generation of entropy is not available, the entropy 
convection equation, as presented in Eq. 5, is useless for computation purposes in this regard. 

Nevertheless, using the entropy field as a diagnosis tool reveals an interesting insight as to whether a solution 
behaves properly or not. A case of violating the second law of thermodynamics in a numerical solution rarely occurs, 
except for the scheme known to branch out to a physically invalid solution, such as expansion shock wave. For such 
a scheme, a fix is usually implemented in a code by adding numerical dissipation. It can be said that satisfying the 
entropy-nondecreasing condition in a numerical solution is generally not a problem; on the other hand, preventing 
the entropy from increasing where it should not is a challenge numerically. One example is the blunt-body carbuncle 
solution by the Roe scheme (with and without entropy fix), where the entropy will be generated, but irregularly 
corresponding to the carbuncle shock profile. Another is the receding flow. In the case of the Roe-flux solution of 
the 'blunt body, knowing that the carbuncle solution does not violate the entropy condition is comforting, but it does 
not help resolve the problem. The entropy condition for this cases is merely a necessary condition. Hence, we shall 
primarily use the entropy function as an additional indicator to evaluate the performance of a numerical flux function. 

IV. A new scheme: AUSM + -ups 

From the above results and discussion, we can arrive at a conclusion that (1) the entropy, hence overheating, is 
generated in the first instance of time evolution in computation; it stays, grows, and reaches a saturated value. The 
problem should be dealt with at the fundamental level-the current discretization framework is incapable of coping 
with it. More precisely, averaging of all waves consisting of nonuniform flow structures, see Fig. 10, is performed at 
the end of each time evolution step, in spite of an exact treatment of flux function as in the case of employing the exact 
Riemann solver. This ultimately leads to the creation of entropy. 


t 

t 
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This also suggests that the spatial averaging of all wave structures to provide the cell-averaged value should be 
avoided, instead a precise point value should be sought. One way to achieve this is through the use of characteristic 
relationships. The decisive drawback of the conventional approach to implementing the method of characteristics is 
based on a floating grid system (namely grid locations change continually with the flow) in order to avoid interpolation 
of flow variables. This obviously makes coding for 3D flows arduously complicated, resulting in the fate of the 
method of characteristics. Another fundamental deficiency is that the characteristic formulation assumes a smooth 
flow where isentropic condition is the inherent limitation. 0 However, this is precisely what we desire for treating 
the overheating problem. The immediate question is how we deal with the grid issue; the conventional approaches 
is clearly not desirable for the reason mentioned above. Instead, we propose to employ the fixed grid system and 
choose to interpolate variables, not coordinates! Then the interpolation becomes easy in this Eulerian grid system. 
Consequently, this framework is simple and general for extension to 3D grids. 

IV.A. Characteristic equations 

The characteristic equations derived from the ID Euler equations of ideal gas are expressed in terms of the eigenvalues 
(wave speeds) A^ = u — a, A^ = w, A^ 3 ^ = u + a: 


dp — padu = 0, 

along 

rH 

II 

R |^ 

(9) 

dp + padu = 0, 

along 

dx ( 3 ) 

dt 

(10) 

ds = 0, 

along 

dx __ ^(2) 
dt 

(ID 


Notice that the last equation is the isentropic condition, similar to Eq. 5, it is a statement stronger than the usual 
one dp/dp = a 2 because the former implies the latter, but the reverse is not necessarily true by incurring discretization 
errors. This equations system allows us to update the flow solution q = (p, p, u ) at £ n+1 from q 71 . A simple upwind 
discretization connects the q at (xi, t n+ 1 ) to q at t n of three distinct points on the corresponding characteristic curves, 
respectively denoted as x^\ x^ 2 \ as shown in Fig. 13. Since they lie on the characteristic curves they are referred 
to as characteristic nodes hereafter. These points generally do not coincide with the grid points, to avoid confusion 
with grid indices we use superscripts to denote the characteristic nodes. A simple discretization produces 


p" +1 — pi 1 ! — cr! 1 l(u” +1 — ul 1 !) = 0, 

along 

dx 

dt 

= AW; 

(12) 

p n+l_ p (3) +(T (3) (u n + l_ u (3) ) = 0j 

along 

dx 

dt 

= A( 3 >; 

03) 

s r +1 =s (2 \ 

along 

dx 

dt 

= A (2) ; 

(14) 


where cr^ = (pa)^ l \ cr^ = (a 2 )( 2 \ cr^ = (pa)^ . The values of q at these points on the characteristic curves are 
determined by a spatial interpolation of q at their neighboring grid points. The key step in the interpolation is to locate 
the characteristic nodes-a straightforward method is proposed in the following. Within the order of accuracy consistent 
with the discretization carried out so far, the characteristic curves are piecewise linear in At = t n+1 — t n . The 
displacement of a characteristic line can be easily determined. Considering computational cell j and its neighboring 
cells, and assuming A^ , A^2 1 > 0, i = 1, 2, 3 for illustration, see Fig. 1 1 . 


\^At = a {i) Ax, cv (i) > 0, i = 1,2,3. (15) 

In other words, the characteristic nodes are defined via the parameters a^ l \ expressing the distance from the grid point 
Xj as a fraction of the cell size. It is known once A^ is determined, as given by 


a£ } = tip - S x xfa^Ax, i = 1, 2, 3 


( 16 ) 


c The nonconservative form of characteristics is also a matter of concern for handling shocks. 


8 of 13 


American Institute of Aeronautics and Astronautics 




|< — Ax — >| 


Figure 11. characteristic curves and nodes. 


Hence, we get the eigenvalues/wave speeds at the characteristic nodes and their locations: 


a2> = 


k«> 


1 + <5 X A^ At 


a (0 = A<*> — 
a Ax’ 


where the gradients J X A^ , i = 1,2,3 can be estimated by either upwind or central differencing. 


A (i) 

r v(i) _ A 3 ~ A 3 - 1 

“ Ax 


or 


\(0 \(») 

X \(0 _ A i +1 ~ A J -1 

^ “ 2Ax • 


(17) 

(18) 


(19) 

( 20 ) 


We have found no clear influences on results to warrant one formula over the other. Then, a linear interpolation of q 
or any variable can be made of (q^-i , q*), e.g.. 


= (1 - a^)pj + oP^Pj- i, 

for k = 1,3; 

(21) 

= (1 — a^)uj -b ot^Uj- 1, 

for k = 1,3; 

(22) 

f ( 2 ) = (1 — ot^)sj -b a^Sj- 1. 


(23) 


Expression for A^ , A^ x < 0 or A^A^ < 0 can be derived in a similar manner. 
Now solving for q n+1 from Eq. 12-14 yields 


„n+l 

Pj 


„n+l 


o n +l 

s j 


1 

(T^ 1 ) ■+■ <7 ( 3 ) 

1 

(jC 1 ) + cr( 3 ) 


(cr^p^ + a^p^ + — i/ 3 ^)), 

+ (jpW — p^)), 


s (2) , giving pj +1 = ( 


P? +1 


exp (s" +1 ) 


>*• 


(24) 

(25) 

(26) 


The receding problem is now solved entirely by the method of characteristics given above and the solution is shown 
in Fig. 13. The overheating is gone! Here we have results of conserving the entropy and eradicating the overheating 
simultaneously; it confirms the validity of the characteristics method and the entropy’s connection to overheating. 
Moreover, the framework for employing the characteristics in a fixed Eulerian grid is simple. 

However, as stated earlier the method of characteristics does not satisfy the jump conditions across a shock wave 
because of its smoothness requirement, a computed result for the Sod problem is shown in Fig. 13; the entropy is 
continuous across the shock wave and as a result the jump condition is not satisfied-higher density and lower pressure 
than the correct values behind the shock. 
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(a) Receding flow. 


(b) Shock tube flow. 


Figure 12. Solutions by the method of characteristics: (a) receding flow and (b) shock tube problem. 


IV.B. AUSM+ups 


Clearly, the characteristics method is not general enough to be used for flows involving shock waves; however it has 
attractive properties for smooth flows. Hence, a proper combination of the conventional Eulerian CFD method with 
the characteristics method will potentially lead to preservation of the best features of both methods. The key lies in 
the statement that facilitates the transition from one to the other. We suggest the following condition for invoking the 
characteristics calculation when 


du 



0 . 


(27) 


This condition can be easily recognized for an expanding flow, such as the receding flow. This condition works for 
several standard cases tested, the findings are all consistent with what was shown above. 
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(a) Receding flow. 

Figure 13. Solutions by the method 


0(h), AUSM+-ups 



(b) Shock tube flow. 

-istics: (a) receding flow and (b) shock tube problem. 


Hence the new method is referred to as AUSM + -ups, as it includes the step for conserving the entropy function 
when appropriate. 
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IV. C. Further verification 

In the following additional tests on problems emphasizing different flow features are included for comparison of 
performance of three methods: AUSM + -up, charateristics, and AUSM+-ups. All calculations presented below are 
carried out on a domain of 200 cells using the first-order explicit Euler method for time integration and CFL=0.5. 

1. Slowly moving contact discontinuity 

The initial condition is set as: 

{p,u,p)l = (0.1, 0. lax,, 1.0), (p,u,p) R = (1.0, 0-la^, 1.0). (28) 

We see in Fig. 14 that all three methods perform rather similarly, all badly smear the contact discontinuity as usually 
observed; for this case the results by the characteristics and AUSM+-upsare identical since the condition Eq. 27 is 
met during the entire computation. 
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Figure 14. Solutions for slowly moving contact discontinuity: (a) AUSM+-up, (b) characteristics method and (c) AUSM + -ups. 


2. Asymmetric receding flow 

We now consider a more complicated receding flow with asymmetric initial conditions between the left and right 
states, which give rise to a contact discontinuity of finite strength is generated. 

(p,u,p) L = (1.0,-2.5a L ,2.0), (p,u,p) R = (1.0, -2.5a L ,0.5). (29) 

The AUSM+-up result, shown in Fig. 15 again shows an entropy generation skewed behind the contact discontinuity, 
thus producing overheating. The results by the characteristics and the new AUSM + -ups are mostly similar-with no 
entropy anomaly and a correct behavior of falling temperature; however three are quantitative differences in the small 
region immediately behind the contact discontinuity. Notice also the smearing at the contact discontinuity is vividly 
manifested by the entropy profile. The pressure and density profiles seem to be in good agreement with the exact 
solution, while the close relationship between the temperature and entropy variations is evident in this problem. 


3. Inverse shock wave 

Now we consider the inverse shock wave case in which the thermodynamic relationship between the left and right 
states satisfy the jump condition corresponding to a given Mach number, except the flow now is moving away from 
the discontinuity, instead of going towards it. The Roe method is known to break down, unless an entropy fix term is 
added. The initial conditions for the right state are: 

(p, u,p)* = (1/6, 5o L , 1/100). (30) 

and the shock relation for M R = u R /a,L = 5 sets (p,u,p)z, = (0.833333,0.0,0.415228). 

Again, we see in Fig. 16 the AUSM + -up result has an entropy generation accompanied with a temperature rise 
(albeit slight but recognizable); the other two solutions are quite similar to each other, with correct behavior in entropy 
and temperature, also smeared across the contact discontinuity as in the case of the above asymmetric receding flow. 
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Figure 15. Solutions for an asymmetric receding flow: (a) AUSM+-up, (b) characteristics method and (c) AUSM + -ups. 
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Figure 16. Solutions for inverse shock problem: (a) AUSM + -up, (b) characteristics method and (c) AUSM+-ups. 


4. Strong colliding flow with a nonvanishing contact 

Finally we test the new method on a shock-contact-shock problem, suggested as problem #5 in Toro’s book 3 (p. 
129), Initially two supersonic streams move towards each other at M=l. 88926 with different sets of thermodynamic 
states, as given below. 

(, p,u,p) L = (5.99924,19.5975,460.894), (p,u,p) R = (5.99924,-6.19633,46.0950). (31) 

The contact discontinuity sandwiched by the two rightward-moving shocks is quite strong; the trailing shock wave 
moves very slowly in comparison with the leading one. 
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Figure 17. Solutions for a shock-contact-shock problem (#5 in Toro’s test problems 3 ): (a) AUSM + -up, (b) characteristics method and (c) AUSM + -ups. 

The characteristics method failed to evolve the pattern and instead maintained the same initial state, its solution is 
not included. The both results by AUSM + -up and AUSM + -ups are very nearly identical, as displayed in Fig. 17. 
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V. Conclusion 


We have been interested in studying some problems that remain intractable with the current state of the art numer- 
ical methods; the overheating phenomenon seen during rarefaction (e.g. in expanding flow) has proven to be tough 
to crack. Some important findings from our previous studies include: (1) the overheating problem remains elusive 
for the current state of the art upwind methods, (2) the maximum error in temperature rise cannot be simply reduced 
by refining the cell or time step size, (3) the limiter function and the set of variables chosen for interpolation affect 
the evolution of the temperature field, (4) no convergence in temperature in the sense of norm (maximum error) 
can be established as the temperature reaches a saturated value, and (5) there is a close correspondence between the 
behaviors of temperature and entropy. 

The objectives of this study are: (1) to further explore and determine the cause of the overheating in the receding 
problem and (2) to ultimately develop a numerical method that is capable of eliminating the overheating problem. We 
have identified the one-to-one link of this temperature anomaly to the entropy generation, which occurs within the first 
time step and grows quickly to a saturated value. A theoretical analysis on the relationship of entropy and temperature 
is developed and validated by the numerical result. It is then concluded that preservation of entropy constancy for a 
smooth flow such as the receding flow is necessary to completely eradicate the overheating problem. The method of 
characteristics is the one that fits this need, but the conventional approach for implementing it is arduously complicated 
because of its continually adjusting the grid. 

In this paper, we propose a new formulation that simply avoids the floating grid concept and is based on a fixed 
(Eulerian) grid, thereby greatly simplifying the implementation of the characteristics method, making it easily extend- 
able to 3D problems. This characteristics method alone is shown capable of maintaining isentropic condition, thus 
completely removing overheating. 

However, the characteristics method is poorly equipped to cope with entropy-increasing discontinuities (shock 
wave in particular), limiting its usefulness. A new method is then developed by hybridizing the conventional finite 
volume upwind method, here the AUSM + -up method, and the just developed characteristics method, with a one-line 
and general logic. The result is that we have a unified approach that can deliver a smooth solution without afflicting 
with an entropy generation and meanwhile can provide a discontinuous solution with correct jumps. The new method 
is thus called AUSM + -ups, as the newest member of the AUSM family. 
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